# Temas para os gráficos
theme.base <- theme_minimal(base_size = 11) +
theme(
axis.text = element_text(size = 8),
plot.title = element_text(hjust = 0.5, size = 11, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, size = 9),
plot.caption = element_text(size = 8),
axis.title = element_text(size = 8),
legend.title = element_text(size = 8),
panel.grid.major = element_line(colour = "grey90", linewidth = 0.5),
panel.grid.minor = element_line(
colour = adjustcolor("grey90", alpha.f = 0.4),
linewidth = 0.5
),
panel.border = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
axis.line.x = element_line(colour = "grey"),
axis.line.y = element_line(colour = "grey"),
)
plotly.base <- function(p) {
p %>%
layout(margin = list(b = 60, t = 80)) %>%
config(mathjax = 'cdn')
}
my.plotly <- function(p) {
if (knitr::is_latex_output()) return(p)
ggplotly(p) %>% plotly.base
}Durante o curso de CEQ, vimos, predominantemente, cartas de controle de Shewhart. Este, no entanto, usa estatísticas apenas da última amostra e ignoram a informação de amostras anteriores, o que a torna pouco sensível para detectar mudanças pequenas, \(\leq1.5\sigma\), no processo.
Para contornar essa limitação, existem duas alternativas principais: CUSUM e EWMA. Neste texto, vamos abordar a segunda.
EWMA é a sigla para Exponential Weighted Moving Average, ou Média Móvel Exponencialmente Ponderada. A ideia é simples: ao invés de usar apenas a última amostra, vamos usar uma média ponderada de todas as amostras anteriores. A ponderação é exponencial, o que significa que amostras mais antigas têm menos peso que amostras mais recentes.
Matematicamente, a média EWMA é dada por:
\[ z_i = \lambda x_i + (1-\lambda)z_{i-1} \]
onde \(z_t\) é a média EWMA no instante \(t\), \(x_t\) é a amostra no instante \(t\) e \(\lambda \in (0, 1]\) é o fator de suavização aplicado a partir da amostra \(i=1\), sendo \(z_0 = \mu_0\).
Desta forma temos que, para o \(\lambda\), quanto mais próximo de 1, mais peso é dado à amostra mais recente. Por outro lado, quanto mais próximo de 0, mais peso é dado às amostras antigas.
Uma vez que o EWMA é uma média ponderada de todas as amostras anteriores, ele é pouco sensível à suposição de normalidade dos dados. Isso o torna ideal para monitorar observações individuais.
Se as observações \(x_i\) são variáveis aleatórias independentes, com variância \(\sigma^2\), então a variância da média EWMA é dada por:
\[ \sigma^2_{z_i} = \sigma^2 \left( \frac{\lambda}{2-\lambda} \right) \left[ 1 - \left( 1-\lambda \right)^{2i} \right] \] Desta forma, os limites de controle para o EWMA são dados por:
\[ LC = \mu_0 \pm L\sigma \sqrt{\frac{\lambda}{2-\lambda} \left[ 1 - \left( 1-\lambda \right)^{2i} \right]} \]
onde \(L\) é o fator de multiplicação dos limites de controle, \(\mu_0\) é a média inicial, \(\sigma\) é o desvio padrão, \(\lambda\) é o fator de suavização e \(i\) é o número da amostra.
Vamos simular um processo com média \(\mu = 100\) e desvio padrão \(\sigma = 5\). Vamos monitorar o processo com um fator de suavização \(\lambda = 0.2\).
ewma <- function(x, lambda, mu, sigma, L) {
ewma <- numeric(length(x))
ewma[1] <- mu
LCS <- numeric(length(x))
LCI <- numeric(length(x))
for (i in 1:length(x)) {
ewma[i] <- lambda * x[i] + (1 - lambda) * ifelse(i == 1, mu, ewma[i - 1])
amostra_sd = sigma * sqrt((lambda / (2 - lambda)) * (1 - (1 - lambda) ^ (2 * i)))
LCS[i] <- mu + L * amostra_sd
LCI[i] <- mu - L * amostra_sd
}
ewma <- data.frame(
observacao = 1:length(x),
medida = x,
ewma = round(ewma, 3),
LCS = round(LCS, 3),
LCI = round(LCI, 3),
fora_de_controle = ewma < LCI | ewma > LCS
)
return(ewma)
}exemplo_medidas_media <- 10
exemplo_medidas_sd <- 1
exemplo_lambda <- 0.1
exemplo_limite_de_controle <- 2.7
exemplo_medidas <- c(
9.45,
7.99,
9.29,
11.66,
12.16,
10.18,
8.04,
11.46,
9.20,
10.34,
9.03,
11.47,
10.51,
9.40,
10.08,
9.37,
10.62,
10.31,
8.52,
10.84,
10.90,
9.33,
12.29,
11.50,
10.60,
11.08,
10.38,
11.62,
11.31,
10.52
)
ewma_exemplo <- ewma(exemplo_medidas, exemplo_lambda, exemplo_medidas_media, exemplo_medidas_sd, exemplo_limite_de_controle)
ewma_exemplomy.plotly(
ewma_exemplo %>%
ggplot() +
geom_line(aes(
x = observacao, y = ewma, color = "EWMA"
)) +
geom_line(aes(
x = observacao, y = LCS, color = "Limite Inferior"
), linetype = "dashed") +
geom_line(aes(
x = observacao, y = LCI, color = "Limite Superior"
), linetype = "dashed") +
geom_hline(
aes(yintercept = exemplo_medidas_media, color = "Média do processo"),
linetype = "dotted"
) +
geom_point(aes(
x = observacao, y = medida, color = "Medida"
), shape = 1) +
geom_point(
data = . %>% filter(fora_de_controle),
aes(x = observacao, y = ewma, color = "Fora de controle"),
size = 2,
shape = 4
) +
labs(
x = "Amostra",
y = "EWMA",
title = "Monitoramento do processo com EWMA",
color = "Legenda"
) +
scale_color_manual(
values = c(
"EWMA" = "blue",
"Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
"Limite Superior" = adjustcolor("red", alpha.f = 0.5),
"Média do processo" = adjustcolor("black", alpha.f = 0.5),
"Medida" = adjustcolor("black", alpha.f = 0.6),
"Fora de controle" = adjustcolor("#f42f2f", alpha.f = 0.8)
)
) +
theme.base
)O \(\lambda\) é um parâmetro importante no EWMA. Ele controla a sensibilidade do gráfico. Quanto mais próximo de 1, mais sensível o gráfico será a mudanças no processo. Por outro lado, quanto mais próximo de 0, menos sensível o gráfico será. Além disso, quando \(\lambda = 1\), teremos a carta de controle de Shewhart, pois a média EWMA será igual à média das amostras.
ewma_exemplo2 <- bind_rows(
ewma_exemplo %>% select(observacao, ewma) %>% mutate(lambda = as.factor(exemplo_lambda)),
ewma(
exemplo_medidas,
0.01,
exemplo_medidas_media,
exemplo_medidas_sd,
exemplo_limite_de_controle
) %>% select(observacao, ewma) %>% mutate(lambda = as.factor(0.01)),
ewma(
exemplo_medidas,
0.5,
exemplo_medidas_media,
exemplo_medidas_sd,
exemplo_limite_de_controle
) %>% select(observacao, ewma) %>% mutate(lambda = as.factor(0.5))
)
my.plotly(
ewma_exemplo2 %>%
ggplot() +
geom_line(aes(
x = observacao,
y = ewma,
color = "EWMA",
linetype = lambda
)) +
geom_line(
data = ewma_exemplo,
aes(x = observacao, y = LCS, color = "Limite Inferior"),
linetype = "dashed"
) +
geom_line(
data = ewma_exemplo,
aes(x = observacao, y = LCI, color = "Limite Superior"),
linetype = "dashed"
) +
geom_point(
data = ewma_exemplo,
aes(x = observacao, y = medida, color = "Medida"),
shape = 1
) +
labs(
x = "Amostra",
y = "EWMA",
title = "EWMA para diferentes valores de λ",
color = "Legenda"
) +
scale_color_manual(
values = c(
"EWMA" = "blue",
"Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
"Limite Superior" = adjustcolor("red", alpha.f = 0.5),
"Medida" = adjustcolor("black", alpha.f = 0.6)
)
) +
theme.base
)Para Montgomery, em geral, valores de \(\lambda\) entre 0.05 e 0.25 são recomendados. Valores menores que 0.05 são muito insensíveis a mudanças no processo, enquanto valores maiores que 0.25 são muito sensíveis a variações normais do processo.
lambdas = expand.grid(lambda = c(0.1, 0.2, 0.3, 0.5, 0.6), amostra = c(0:10)) %>%
mutate(
peso = round(lambda * (1 - lambda) ^ amostra, 4),
lambda = as.factor(lambda)
)
my.plotly(
lambdas %>%
ggplot(aes(
x = amostra, y = peso, color = lambda
)) +
geom_line() +
labs(
x = "Idade da amostra",
y = "Peso",
title = "Peso das amostras em função do fator de suavização",
color = "Valor de λ"
) +
theme.base
)O EWMA também pode ser usado para monitorar proporções. Neste caso, a estatística EWMA permanece a mesma:
\[ z_i = \lambda x_i + (1-\lambda)z_{i-1} \]
onde, agora, \(x_i \sim Poi(l)\) é a contagem na amostra \(i\), com \(z_0 = \mu_0\) a taxa em controle.
Já o limite de controle é dado por:
\[ \begin{aligned} \text{LCS} &= \mu_0 + A_S \sqrt{\frac{\lambda\mu_0}{2-\lambda} \left[ 1 - \left( 1-\lambda \right)^{2i} \right]} \\ \text{LCI} &= \mu_0 - A_I \sqrt{\frac{\lambda\mu_0}{2-\lambda} \left[ 1 - \left( 1-\lambda \right)^{2i} \right]} \end{aligned} \]
onde \(A_S\) e \(A_I\) são os fatores de multiplicação para o limite superior e inferior, respectivamente. Muitas vezes, \(A_S = A_I = A\).
ewma_atributos <- function(x, lambda, mu, A) {
ewma <- numeric(length(x) + 1)
ewma[1] <- mu
LCS <- numeric(length(x))
LCI <- numeric(length(x))
for (i in 1:length(x)) {
ewma[i + 1] <- lambda * x[i] + (1 - lambda) * ewma[i]
amostra_sd = A * sqrt((lambda * mu) / (2 - lambda) * (1 - (1 - lambda) ^ (2 * i)))
LCS[i] <- mu + amostra_sd
LCI[i] <- mu - amostra_sd
}
ewma <- data.frame(
observacao = 1:length(x),
medida = x,
ewma = round(ewma[-1], 3),
LCS = round(LCS, 3),
LCI = round(LCI, 3),
fora_de_controle = ewma[-1] < LCI | ewma[-1] > LCS
)
return(ewma)
}set.seed(4)
exemplo_medidas_mu <- 0.1
exemplo_medidas_atributos <- unlist(replicate(30, rbinom(1, 1, 0.1), simplify = FALSE))
exemplo_lambda_atributos <- 0.1
exemplo_limite_de_controle_atributos <- 2.7
ewma_atributos_exemplo <- ewma_atributos(
exemplo_medidas_atributos,
exemplo_lambda_atributos,
exemplo_medidas_mu,
exemplo_limite_de_controle_atributos
)my.plotly(
ewma_atributos_exemplo %>%
ggplot() +
geom_line(aes(
x = observacao, y = ewma, color = "EWMA"
)) +
geom_line(aes(
x = observacao, y = LCS, color = "Limite Inferior"
), linetype = "dashed") +
geom_line(aes(
x = observacao, y = LCI, color = "Limite Superior"
), linetype = "dashed") +
geom_hline(
aes(yintercept = exemplo_medidas_mu, color = "Média do processo"),
linetype = "dotted"
) +
geom_point(aes(
x = observacao, y = medida, color = "Medida"
), shape = 1) +
geom_point(
data = . %>% filter(fora_de_controle),
aes(x = observacao, y = ewma, color = "Fora de controle"),
size = 2,
shape = 4
) +
labs(
x = "Amostra",
y = "EWMA",
title = "Monitoramento do processo com EWMA para atributos",
color = "Legenda"
) +
scale_color_manual(
values = c(
"EWMA" = "blue",
"Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
"Limite Superior" = adjustcolor("red", alpha.f = 0.5),
"Média do processo" = adjustcolor("black", alpha.f = 0.5),
"Medida" = adjustcolor("black", alpha.f = 0.6),
"Fora de controle" = adjustcolor("#f42f2f", alpha.f = 0.8)
)
) +
theme.base
)Além de monitorar processos, o EWMA também pode ser usado para prever valores futuros. Desta forma, pode ser usada como base de um processo de controle dinâmico.
Ou seja, a média EWMA pode ser usada como preditor para o próximo valor do processo, sinalizando quando o processo irá sair de controle. Além disso, a diferença entre o valor observado e o valor objetivo pode ser usada para determinar o quanto deve ser ajustado.
Assim, o valor predito é dado por:
\[ \hat{z}_{i} = z_{i-1} + \lambda_1 e_i + \lambda_2 \sum_{j=1}^{i} e_j + \lambda_3 \nabla e_i \]
onde, \(e_i = x_i - z_{i-1}\) é o erro na previsão, \(\nabla e_i = e_i - e_{i-1}\) é a primeira diferença entre os erros e \(\lambda_1\), \(\lambda_2\) e \(\lambda_3\) são os fatores de ponderação escolhidos tais que resultam na melhor performance do preditor.
exemplo_medidas_media <- 10
exemplo_medidas_sd <- 1
exemplo_lambda <- 0.1
exemplo_limite_de_controle <- 2.7
lambda1 <- 0.04 # peso do erro de previsão (proporcional)
lambda2 <- 0.01 # peso da soma cumulativa dos erros (integral)
lambda3 <- 0.10 # peso da diferença de erros (diferencial)
medidas_valores <- exemplo_medidas[1:(length(exemplo_medidas) - 2)]
ewma_valores <- ewma(medidas_valores, exemplo_lambda, exemplo_medidas_media, exemplo_medidas_sd, exemplo_limite_de_controle)
# prever o próximo valor
predicoes <- numeric(length(medidas_valores))
predicoes[1] <- exemplo_medidas_media
predicoes[2] <- ewma_valores$ewma[2]
erro_previsao <- 0
diff_erros <- 0
soma_erros <- 0
for (i in 3:length(medidas_valores)) {
erro_previsao <- medidas_valores[i] - predicoes[i - 1] # Calcula o erro atual
diff_erros <- erro_previsao - (medidas_valores[i - 1] - predicoes[i - 2]) # Diferença de erros
soma_erros <- soma_erros + erro_previsao # Atualiza a soma cumulativa dos erros
predicoes[i] <- predicoes[i - 1] +
lambda1 * erro_previsao +
lambda2 * soma_erros +
lambda3 * diff_erros # Calcula a previsão
}
predicoes_ = c(NA, predicoes)
ewma_preditor <- bind_rows(
ewma_valores,
data.frame(
observacao = ewma_valores$observacao[length(ewma_valores$observacao)] + 1,
medida = NA,
ewma = NA,
LCS = ewma_valores$LCS[length(ewma_valores$LCS)],
LCI = ewma_valores$LCI[length(ewma_valores$LCI)],
fora_de_controle = NA
)
)
ewma_preditor$predicao <- predicoes_
ewma_preditor$fora_de_controle_predicao <- (
ewma_preditor$predicao < ewma_preditor$LCI | ewma_preditor$predicao > ewma_preditor$LCS
)my.plotly(
ewma_preditor %>%
ggplot() +
geom_line(aes(
x = observacao, y = ewma, color = "EWMA"
)) +
geom_line(aes(
x = observacao, y = predicao, color = "Predição"
), linetype = "dotted") +
geom_line(aes(
x = observacao, y = LCS, color = "Limite Inferior"
), linetype = "dashed") +
geom_line(aes(
x = observacao, y = LCI, color = "Limite Superior"
), linetype = "dashed") +
geom_hline(
aes(yintercept = exemplo_medidas_media, color = "Média do processo"),
linetype = "dotted"
) +
geom_point(aes(
x = observacao, y = medida, color = "Medida"
), shape = 1) +
geom_point(
data = . %>% filter(fora_de_controle),
aes(x = observacao, y = ewma, color = "Fora de controle"),
size = 2,
shape = 4
) +
geom_point(
data = . %>% filter(fora_de_controle_predicao),
aes(x = observacao, y = predicao, color = "Fora de controle (predição)"),
size = 2,
shape = 4
) +
labs(
x = "Amostra",
y = "EWMA",
title = "Monitoramento do processo com EWMA e predição",
color = "Legenda"
) +
scale_color_manual(
values = c(
"EWMA" = "blue",
"Predição" = "darkgreen",
"Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
"Limite Superior" = adjustcolor("red", alpha.f = 0.5),
"Média do processo" = adjustcolor("black", alpha.f = 0.5),
"Medida" = adjustcolor("black", alpha.f = 0.6),
"Fora de controle" = adjustcolor("#f42f2f", alpha.f = 0.8),
"Fora de controle (predição)" = adjustcolor("#f42f2f", alpha.f = 0.8)
)
) +
theme.base
)De uma forma geral, o CUSUM possui mais poder que o EWMA para detectar mudanças pequenas no processo. No entanto, o EWMA é mais simples de implementar e interpretar.
No gráfico a seguir é mostrado o resultado da otimização por Monte Carlo dos parâmetros CUSUM e EWMA. Nota-se que o CUSUM é mais sensível às mudanças no processo.
comparacao_ewma_cusum <- readRDS("estatisticas-combinadas-resumo.rds")
pp_ <- my.plotly(
comparacao_ewma_cusum %>%
ggplot() +
geom_line(aes(
x = h1_phi,
y = mean,
color = parametro,
linetype = algoritmo
)) +
geom_point(aes(
x = h1_phi,
y = mean,
color = parametro,
shape = algoritmo
)) +
geom_vline(
xintercept = 0.2,
linetype = "dotted",
color = "gray70"
) +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
linetype = "Algoritmo",
title = "Fração de pontos fora de controle",
color = "Valores dos parâmetros",
shape = "Algoritmo"
) +
theme.base
)
if (!knitr::is_latex_output()) {
pp_ %>% layout(annotations = list(
x = 0.2,
y = 0.2,
text = "Processo sob controle",
textangle = -90
))
} else {
pp_
}Apesar disso, o ewma possui uma variabilidade menor que o CUSUM.
comparacao_ewma_cusum <- readRDS("estatisticas-combinadas.rds")
pp_ <- my.plotly(
comparacao_ewma_cusum %>%
ggplot() +
geom_boxplot(
aes(
x = h1_phi,
y = fracao_fora_de_controle,
fill = parametro,
shape = algoritmo
),
) +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
fill = "Valores dos parâmetros",
linetype = "Algoritmo",
title = "Fração de pontos fora de controle",
color = "Valores dos parâmetros"
) +
facet_wrap(~ algoritmo) +
theme.base
)
if (!knitr::is_latex_output()) {
pp_ %>% layout(boxmode = 'group')
} else {
pp_
}